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Abstract 

The dynamic equations of motion of an n-body pendulum with spherical joints are derived to be a mixed system 
of differential and algebraic equations (DAE’s). The DAE’s are kept in implicit form to save arithmetic and 
preserve the sparsity of the system and are solved by the robust implicit integration method. At each solution 
point , the predicted solution is corrected to its exact solution within given tolerance using Newton ’s iterative 
method . For each iteration, a linear system of the form J AX = E has to be solved. The computational cost for 
solving this linear system directly by LU factorization is 0(n 3 ), and it can be reduced significantly by exploring 
the structure of J. This paper shows that by recognizing the recursive patterns and exploiting the sparsity 
of the system the multiplicative and additive computational costs for solving J AX = E are 0(n) and 0(n 2 ), 
respectively. The formulation and solution method for an n-body pendulum is presented. The computational 
cost is shown to be nearly linearly proportional to the number of bodies. 


1 INTRODUCTION 

The general modeling and formulation of an open-chain multi-body system with spherical joints was presented 
by Chou, Singhal, and Kesavan in 1986 [1]. In this paper, we are interested in a single open kinematic chain 
without branching which is a special configuration of a general open-chain system. Much attention was paid to 
this open-chain system by researchers, such as Armstrong [2], because the system is simple and its configuration 
is similar to robot arms. In Armstrong’s work, he presented an O(n) algorithm for the computation of robot 
forward dynamics. However, in the conclusion of his paper he stated that his program worked only for the 
joints with three degrees of freedom (i.e., spherical joints). He also claimed that his algorithm can be enhanced 
to include prismatic and revolute joints. In the author’s opinion, once we have joints which possess less than 
three rotational degrees of freedom connecting two bodies, the bodies are rotationally dependent. In this case, 
the equations of motion can not be decoupled easily, and we may not be able to find an O(n) algorithm for the 
complete calculation of robot forward dynamics unless we use very simple and primitive integration methods 
or other techniques such as parallel computing. In other words, if we use an unreliable integration routine we 
may get invalid solutions. 

A series of n rigid bodies joined sequentially by spherical joints form an n-body pendulum. The bodies are 
allowed to rotate freely in space, but the motion of translation between adjacent bodies is constrained by the 
joint. This sequence of bodies can swing in any direction with one end fixed at the ceiling through a spherical 
joint. Here, we derive the equations of motion of this pendulum in the form of a mixed set of differential and 
algebraic equations (DAE’s) and solve the equations using the robust implicit integration method developed by 
Petzold [3,4]. The DAE’s are kept in implicit form to save arithmetic and preserve the sparsity of the system. 
At each solution point, the predicted solution is corrected to its exact solution within the given tolerance using 
Newton’s iterative method. For each iteration, a linear system of the form J AX = E has to be solved. The 
computational cost for solving this linear system directly by LU factorization is 0(n 3 ). In order to reduce 
computations, we explore the structure of J and take advantage of the system sparsity. This paper shows that 
by recognizing the recursive patterns and exploiting the sparsity of the system the multiplicative and additive 
computational costs for solving J AX = E are O(n) and 0(n 2 ), respectively. 
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2 MATHEMATICAL MODEL 


An n-body pendulum was modeled with graphs by Chou and was presented in [5]. Based on this graph-theoretic 
model, the mathematical model was derived as a set of DAE’s and was written symbolically as 


E(X,X,<) 


where the vector of unknown variables is 


X = [ R* VI Pj Sf ] 3 


( 1 ) 


( 2 ) 


The vectors R& and V& are the collection of displacements and velocities of all the bodies, and the vectors P* 
and S& are the orientation parameters (we use Euler parameters here) and their derivatives. The superscripts 
”u” and ’V’ indicate the dependent set and the independent set, respectively. The equations include the 
following six sets of equations: 

' 'F(V 6 , p 6 , S*, S», t) 

*g(r*. p»; t) 


E = 


T I 2 (P*. S h , t) 

[ 'I 3 (PJ, SJ, t) 

and they are a total of 14n scalar equations in 14n unknown variables. 


‘G 

"T 


(Vs, Pt, 

(P», 


Ss, tj 


= 0 


(3) 


3 SOLVING THE FULL SYSTEM 

At each solution point, we can solve the system equations E directly, using the implicit integration method 
[3,4] where the predicted solution is corrected to its exact solution within a given tolerance using Newton’s 
iterative method. For each iteration, a liner system 


J AX = E 


has to be solved by LU factorization. The matrix J is the system Jacobian matrix which is specified by the 
formula given by Petzold [4]. The vectors AX and E are the vectors of corrections and residuals, respectively. 
Letting the Jacobian of an implicit function F with respect to the variable V be Fy , we can write the Jacobian 
matrix J and the vectors AX and E of the full system symbolically as follows: 


l G R 

0 

0 

0 

0 


0 

0 

0 

"Ti 


r r2 
l S u 

0 

0 


0 

*G P . 


' AR» 


‘G I 

*<3 S . *G P u 

*G P , 


AV* 


‘G 

0 0 

r rJ 

J pm 


AS? 


'I 2 

0 

rjS 

1pm 


ASJ 

— 

'I 3 

0 r I l p . 

r rl 

1pm 


AP? 
A ZA 


'I 1 

r F s * r Fp. 

*Fp. 


. AP S . 


T J 


(4) 


Solving the full system means that the system is solved by implicit integration in which the full linear system, 
J AX = E, is solved by LU factorization without reducing its size. However, some structural information in J 
and E can be utilized to reduce computations when they are evaluated. 


3*1 Recursiveness and Special Structures 

Some equations in E possess recursive properties which can be further utilized to reduce computational cost. 
This recursive nature also causes the sub- matrices in J to have special structures which can be exploited to 
minimize the computational cost. 


3.1.1 Translational Kinematic Constraints *G 
In [5], the translational kinematic equations were derived as 

i 

*Gi = R». — - Aiii + ^ Ak&k = 0 (5) 

k=l 
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Let 




— T!i-i Ak&k 

Xi = Ai&i 

Ao = 0 


then, (5) becomes *G< = R &i — R jpi — + A< = 0, or recursively 

*G< = Rjj — R, Pl — A{Ti + (Aj.i+ii) = 0 


( 6 ) 


When we evaluate the residual vector *G, we only have to compute R> 4 , R, Pl , A^r, , and Xj for i — 
once. Instead of computing for each equation, we only compute x< and add the previous A,i ; that is, 

A» — A<_i + Xj. 

The Jacobian entries corresponding to the equation *G are 1 Gr and l Gp. The Jacobian of f G correspond- 
ing to R» can be shown to be a unit matrix easily. Here we show that *Gp has a special structure. It is written 

55 T *1^1 **G, d*Gi 

«P», «P»» 

a‘G, a*G» g*G. 

*P», »P», »P»„ 


'G P = 


d‘G 




, «‘G. «‘G. »‘G„ 


( 7 ) 


The matrix in (7) becomes 


l G P = 


and p >4 , 

for each 

5‘G< 

n. L . 

& P*» 

— 0; * 

a*G. 

o 

TpT 

6'G, 

»‘G, 

TpC 


a«G. 
. «p*, 

*«Q. 

®9h 


( 8 ) 


0 

0 


0 

0 

a‘G„ 


( 9 ) 


and it is a blocked lower-triangular matrix. Each blocked entry is a 3 x 4 matrix. The entries in each column 
i (t = !,...,«) in (9) are derived, in Appendix D.l in [5], as 


I 


= 2i Ei a< — 2J5,- ti 
^ = 2Ei a<; k = i + 1 n 


( 10 ) 


Hence, (9) becomes 


'G P = 


2fJ 1 (a 1 - r x ) 0 0 

2E\ 2 JEr 2 * 2 ) 0 

2E\ a x 2 jE?2 *2 2^a(a s - r 3 ) 


2 JE? i ai 


2 JE?2 a 2 


2E$ a 3 


0 
0 
0 

2En(*n ~ in) J 


( 11 ) 


To compute t Gp 1 we only have to evaluate 2Ei a*, 2E{ r*, and 2Ei(*i — r*) once for i = 1, . . . , n. 

The matrix t Gp can be further partitioned into l Gp~ and *Gp. as in (4). This involves permuting selected 
columns; however, the special structure is retained for these matrices after they are partitioned. 
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3.1.2 Translational Velocity Kinematic Constraints *G 
Translational velocity constraints are written as 


n 

‘G< = V ti - AiXi + = 0 


( 12 ) 


k=l 


Let 


{ 


= £*=i^* a * 

i i = Ai*i 
Ao = 0 


then, (12) becomes *G< = V>. - = 0, or recursively 

‘G < = V 4i - -*i»« + &_!+*<) = 0 


(13) 


Similar to (6), we only have to compute Vj,, and &i once. Instead of computing A <, we compute A<_i +i< 
recursively to save computations. 

The Jacobian matrices correspond to ( G are t Gy, 'Gj, and *Gp in which *Gy can be shown to be a unit 
matrix easily, and *Gs and *Gp will be shown to possess special structures. The Jacobian matrix of ‘G with 
respect to S& is defined as 


*G S = 


0‘G 

as h 


»«6. g*G. 
g'G, e-G, 

• ^ 8 k a 


5 *G, 

Ts^ 1 

f>*Ga 

TsT* 


* 8 »i d 8 * a d 8 *» 

Since *6* is a function of s^, s* 3 , . . and s h . } for each row % (t = 1, ...,n) we have 


- = 0; As = * + 1, ...,n 


Hence, (14) becomes 


‘65 = 


da^ 


0 

*«G, *«G, 

0 »i x ”SbTT 


0 

0 


0 

0 


^*G n A *G,» I 
dS hl d S >3 ?S» n J 

The entries in each column t (i = 1, n) in (16) are derived, in Appendix D.l in [5], as 




= — 2 Ei ti 


— 2-®» a<; i = * + !,...,» 


Comparing (17) with (10), we find that 


f G s = 


9*G 


d‘G 


= 'G P 


6 S» d P» 

Computing *<3$ is not required. Once we compute t Gp, we get *Gs. 


(14) 


(15) 


( 16 ) 


(17) 


(18) 
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The Jacobian, t G P , is different. Since ‘G, is a function of p Sl , p >a , .... and p»., *G P is a blocked 
lower- triangular matrix like % Gs and is written as 




d P» 


g*G, 

0 

0 

0 


ap», 

g‘G, 

«‘G, 

0 

0 


»P‘. 

*P», 





o‘G n 

^ . 

a*G n 


*P», 

^P», 

«P*5 

«P*„ - 


(19) are 

derived, 

in Appendix D.l in [5], as 


(19) 


I 


= 2 E t a< - 2 Ei n 
= 2 Eiii\ k = i + 1 » 


( 20 ) 


Hence, (19) becomes 


‘Gp = 


2i7i(ai ~ r x ) 0 0 

2E\ 2J57 2 (a 2 — r 2 ) 0 

2i?i a* 2i? 2 a 2 — * 3 ) 


0 

0 

0 


2Ei &i 


2E 2 a 2 


2E S a 3 


2 i? n (a n - r n) 


( 21 ) 


To compute t Gp } we only have to evaluate 2 Ei a 217* r*, and 217* (a* — r*) once for t — 1, . . . , n. 

The matrices t Gs and l Gp can be further partitioned into l Gs t ^s»> and l Gp« as in (4). This 

involves permuting selected columns. However, after permuting columns the special structures are retained for 
these partitioned matrices. 

3.1.3 Torque-Balance Equations r F 

The torque-balance equations for an n-body pendulum were derived as 


T F < = T’ &i - iiAj Mi(V hi - s) + *iAf ^Af*(V 6k - g) = 0 

k—i 


( 22 ) 


Let 

= E n k=i - g) 

*< = - g) 

r„+i = 0 

then, we can write (22) as r F< = T*. - r {AT + a { Af r { = 0, or recursively 

r Fi = Tj. - iiAfti + &iAj(zi + r i+l ) = 0 

When we compute the residual of r F, we only have to evaluate Tj . , z* , r«.A.f and a < Af r< once for i = 1, n. 

Instead of computing v ,• for each equation, we compute Zj and accumulate it recursively to save computations. 
The recursive term, Zj + t<+i, runs in backward sequence; that is, i = n, ..., 1. 

The Jacobian of 'F has three parts: T Fy, T Fs, and r F P . They also possess special structures. First, the 
Jacobian of 'F with respect to Vs is written as 


(23) 


d'F d T F d T F 

Fy : + <T ; — — <T # — — <T 

dv h dv b dx b 


■ g*F, 

g’F, 

g-F, ' 

»v„ 

«v l3 

*v. n 

g'F, 

g'Fi 
1 — 

.. g'Fi 


oV, 


g'F„ 
1 — “■ 

o 'F n 

.. » 

•v., 


J 


(24) 
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Since r F» is a function of V^, V* i+1 , . . and for each column »(* = !, n) we have 


<r — = 0; + n 

dv bi 


The matrix (24) becomes a blocked upper- triangular matrix and is written as 


T F V = cr 



0 


fl’T, 

a r Fi 



a-F, 

a-F, 



0 

a-F, 




The entries of each row i (i = 1, n) in T Fv are derived as 





= — aMiXiA f + <rMikiAf 
= aMkiiAj; k = i + l, ..., n 


Hence, (26) becomes 


r F v = <r 


Mi aj Af 

0 Ms(& 2 — f 2 )Aj 

0 0 


M 3 inAT ••• 
Aj 

M 3 (o 3 — T 3 )Ag 


0 


0 


0 


M n a. iAT 
M n *tA* 
•Mn&aA 3 

*n)A n m 


( 25 ) 


(26) 


(27) 


(28) 


To compute T Fy, we only have to evaluate <rMj(a,- — i{)Af and aM^iiAf once for i = 1, 

The Jacobian of T F with respect to Pj is similar to (24). However, since r Fj is a function of > the 
olf-diagonal entries become sero: 


d r F k a r F { 


= 0; It = i + 


for i = 1, ..., n. Hence, r Fp becomes a diagonally blocked matrix and is written as 


where 


r F P = 


d T v 


a *F i 

® P‘i 


a-F, 


0 

0 

»P*. 


g r F< 

d P» 4 


d _3i _ -SMsl , r giiLli) 

0 P» 4 ‘ 5 Ps 4 d P»i 

0Tj 

d P hi 


— 2iiGi %i -I- 2a*G t r< 


= {2J < G i + 4GiGjl t Gi - 4(sf iP>i )J < G i - 2^G,} 

— 2fjCr{ Zj + 2&iG{ Tj 

where = 2JiG*p*. and ^ = i/> x are defined in Appendix D.2 in [5]. 


(29) 


(30) 


(31) 
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T Fs has a structure identical to that in r Fp. The ofF-diagonal entries are zero: 


6 T F k d T F k 

+ < 7- 


8 T Fi d r F i _ . . , , 

+ <r— — = 0; k = t + 1,..., n 


d s». d s*. da ih ' d it h 
for i = 1, n. The Jacobian T Fs is diagonally blocked and is written as 

0 


( 32 ) 


r a-F, , s'F. 

as,, ^ e s,, 


r F s = 


(TF (TF 

ds» + a ds h 


a-F, , _aiFi 


+ <r^£» 
a s,„ 


(33) 


where 


^s 6i 


«'F 4 d r T>, a r T &i 

— + <T 


(34) 


d sj, d sj, 5 sj, 

= {4G<Gf I<G< - 4(s£p bi )IiGi + 2&G<} - <r(2J i G i ) 

The derivation of (31) and (34) is given in Appendix D.2 in [5]. 

3.1.4 Equations Concerning Euler Parameters: r I 1 , r I 2 , and r I 3 

When we compute the residuals of r I 1 , r I 2 , and r I 3 , there is no applicable recursive structure that can be 
utilized. However, their Jacobians do have some special structures because these equations are derived for each 
body. Matrices with diagonally blocked structures are to be expected. 

First, the normality constraint for each set of Euler parameters is written as 

(35) 


r n = -i = o 


Since r 7/ is a function of p> i only, we have 


8 r l } 


. = = 0; * = i + l,...,» (36) 

O P b h O p hi 

From (36), the Jacobian matrix of r I 1 with respect to p^, p& 3 , and p* n must be a diagonally blocked 
matrix. The same argument can be applied to r I 2 and r I 3 such that their Jacobians with respect to and 
S b are diagonally blocked. 

Since the blocks concerning these three sets of equations are locally processed, as shown in [5], the Jacobian 
matrices r Jp*, r Jp», and T Jp. possess a similar structure which is illustrated as follows: 

0 

0 €j 


r t T 2 

¥ 


r to 
_Kp^ 

T l\. 



(37) 


where 


[<k e< /«] 

—<r 0 

0 —o' 

0 0 

[ai bi c<] 


0 

0 

—O' 


are defined in Appendix B.l in [5]. 
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4 SOLVING THE REDUCED SYSTEM 


Solving the reduced system means to solve the linear system using the technique of dimension reduction. 
Partitioning the linear system, J AX = E, according to the partition in (4) gives 


Jn 

J 12 

' AXi ‘ 


' Ei ' 

J 21 

J 22 

ax 2 


. E * . 


Since J\\ is a unit upper- triangular matrix and is non-singular, we can obtain 


J 11 J 12 

’ AXi ' 

' Ex ' 

0 Jr 

ax 2 

E R 


For the independent corrections AX 2 , we solve 


by LU factorization where 


The dependent corrections AXi are 


JrAX.2 = Ej* 

(38) 

\ Jr — J 22 — J*\J\\Jvi 

(39) 

\ Er = E 2 — (J2lJii)^l 

obtained by backward substitutions: 


AXx = J 11 (Ei — J 12 AX 2 ) 

(40) 


4.1 Solving for Independent Corrections 

The independent corrections AX 2 are obtained by solving the linear system Jr AX 2 — E r. The evaluation 
of J R and E# is accomplished by developing a set of formulae such that the sparsity of the linear system is 
completely utilized. From (4), we derive the following formulae for Jr and Er: 

Jr = - F P . - ( r .Fv)(‘Gj».) - (*# 4 )('4.) 

-r* 5 )( r is>.) - r*e)ri}>.) ( 41 ) 

and 

E* = 'F - ('fV)‘G - ( f # 4 ri 2 - ( '#$) r I 8 - ( r JPe) (42) 

where 

r - ('f v )(«6,.) 

{ = 'F s . - ('FvK'Gs.) (43) 

l = *F P . - ('FvK'Gp.) 

When we compute Jr and E r using the above formulae, several matrix multiplications will be performed. 
The special structures discussed in the previous section can be exploited to reduce computation cost. 

4.2 Solving for Dependent Corrections 

Once we obtain the independent corrections, AX 2 , the dependent corrections, AXi, can be computed by 
Jn( Ei- J 12 AX 2 ). However, using this formula directly is not practical since Jn has a very simple structure 
with many zeros. In order to utilize the sparsity completely, a set of formulae are derived symbolically for the 
dependent corrections. They are as follows: 


APj = 'i 1 - ( r Ji.,)Ap; 

(44) 

as: = 'i 8 - r4-)AP5 

(45) 

AS“ = 'I 2 - r4.)APJ 

(46) 

AV k = ‘G - (*G S -) AS? - (‘G 5 .)ASJ 


- (*6p.)APr - CGp.)APl 

(47) 

AR t = ‘G - (‘ Gp«) APJ - (*G,.)APJ 

(48) 


Again, when we compute the dependent corrections using the above formulae, the special structures discussed 
in the previous section can be exploited. 
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5 LINEAR COMPUTATIONAL-COST SCHEME 


Rewriting the full system (4) by 

• combining S£ and SJ to form S in the original order 

• combining PJf and PJ to form P in the original order 

• combining r I l and r I 3 to form I 

• combining — 2<x r I 2 (scaled by — 2<r) and *F to form F 

• re-arranging the columns and rows 

• dropping the subscripts and superscripts for clarity and simplicity 


gives 


or 


Gr 

G P 

0 

0 


■ AR 


■ G 

0 

Ip 

Is 

0 


AP 


I 

0 

F P 

F s 

Fy 


AS 


P 

0 

G P 

G s 

Gy 


. AV . 


. G 

r u 

Gp 

0 

0 - 


AR ‘ 


■ G 

0 

U 

Is 

0 


AP 


I 

0 

F P 

F s 

F v 


AS 


P 

. 0 

G P 

G S 

U . 


AV 


. G . 


(49) 


(50) 


where Gr and Gy were shown to be unit matrices. In order to make Ip a unit matrix, Gaussian elimination 
has to be applied to Ip, Is> and I locally. An example for one body is given in Appendix B.2 in [5]. We intend 
to solve this linear system with a computational cost which is linearly proportional to the number of bodies in 
the system. 


5.1 Basic Theory 

Applying the technique of dimension reduction to (50) gives 


$ 

O 

O 


' AR ' 


' G ' 

0 U Is 0 


AP 


I 

0 0 Jp 0 


AS 


E p 

0 0 0 Jr 


AV 


. . 


where 

r J R = U - (65 — Gpls) Jf l F v 
\ Ej, = G - (Gp)I - {Gs-Gpls) Ji'p-Fpl) 


{ J p — F$ — I 9 p Is 

Ej? = F - (F P ) I - (F v ) AV 

Instead of solving (50) directly by LU factorisation, we solve 


JrAV = Ejt 


to obtain AV; then we solve 

I J j?AS == Eji> 


to obtain AS. The rest of corrections can be computed by 


/ AP = I - (I 5 ) AS 
\ AR = G - (G P ) AP 


(51) 

(52) 

(53) 

(54) 

(55) 

(56) 


96 




5.2 Invert ibility 


In [6], we have shown that Jr is non-singular if J is non-singular. From (52), we have to find Jp in order to 
evaluate Jr. We would like to show that Jp is invertible when a — ♦ oo. 

The matrix Jp is diagonally blocked because it is computed by F$ — Fp I 5 where F5, JFp, and 1 5 are 
sill diagonally blocked matrices. Each block at the diagonal position in Jp is computed by its corresponding 
block in JF5, Fp y and 1$. Hence, we write 


Jf* = Fsi ~ FPilSi 


(57) 


for each diagonal block. From (B.12) in Appendix B given in [5], we find that J5. is a matrix scaled by 
Hence, we may write 

F Pi i Si = -r, 

<T 

Also, from (34) we find that Fsi may be written as 

F Si = - <r(2IiGi) + r< 


Hence, we can write 

J Fi = - <r(2 IiGi) + A - -Y { (58) 

c 

If <r — ► 00, we have 

J Fi * - <r(2 IiGi) (59) 

However, we have added one more equation into the F block as mentioned previously. This increases the 
number of equations in each block in F from three to four. The equation is r If scaled by —20-; that is, 


— 2<r T lf = — 2 <rpf s,- = 0 


( 60 ) 


The Jacobian corresponding to s is 

- 2<rpf (61) 

Combining it with (59) gives 

J Fi » -<r(2J,p;) (62) 

Adding equation (60) into block F is equivalent to using Euler’s equations in quaternion space [7]. We 
have proved that the coefficient matrix with respect to s in Euler’s equations for inertial torque in 4^space is 

a non-singular matrix [7]. This matrix is exactly the same as 2Ji(p<) T in (62). Therefore, Jp. as well as Jp 
are invertible if <7 — ♦ 00. 


5.3 Summary of Computational Cost 

The detailed discussion of the linear-cost scheme was presented in [5]. In terms ofmultiplicative(x) and additive 
(+) operations, the computational cost required for solving the reduced linear system for the linear-cost scheme 
is summarized also in detail in [5] as Tables 2 to 6. Adding up the totals in the tables gives the grand total of 
operations required for the linear-cost scheme: 


/ x : 657n - 324 
\ + : 15n 2 + 574n - 321 


(63) 


This cost is not for obtaining one solution point; it is the cost of solving the reduced linear system for one 
Newton’s iteration. The computational cost for reaching a solution point depends on output step size and the 
number of iterations. 


6 IMPLEMENTATION AND RESULTS 

A program, called LINPEN (LINear-cost scheme for simulating an n-body PENdulum), was coded to sim- 
ulate a pendulum. The program is equipped with modules for reporting various physical quantities such as 
displacement, velocity, and acceleration and was run on a VAX 750 computer. The motion of a user-specified 
pendulum can be displayed on a GRINNELL display terminal. 

Each body is drawn as a standardized 3-D polygon and is projected in perspective, on the terminal screen 
according to its simulated position and orientation. An example of simulating a 3-body pendulum for ten 
continuous positions is illustrated in Figure 1. The figure looks exactly the same as displayed on the display 
terminal. 
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6.1 Comparative Study in Computational Cost 

For each solution point, a linear system, J AX = E, may have to be solved by LU factorization several times. 
There are four schemes available to accomplish this: 

• solving the full linear system by LU factorization where J is approximated by numerical difference 

• solving the full linear system by LU factorization where J is evaluated by the exact Jacobian matrix of 
E 

• solving the reduced linear system by LU factorization where J is evaluated by the exact Jacobian and 
is further reduced to its optimal size using the technique of dimension reduction 

• solving the reduced linear system by LU factorization where J is evaluated by the exact Jacobian and 
is further reduced to its optimal size using the linear-cost scheme 

For each iteration, the analytical count of operations required to solve the full system, reduced system, and 
the linear-cost system with exact Jacobian is tabulated in Table 1. Although the additive operations for the 
linear-cost system is 0(n 2 ), we may consider the cost linear because the amount of time required to perform a 
multiplication or division on a computer is about the same and is considerably greater than that required to 
perform an addition or subtraction. 

The simulation of an n-body pendulum was run for the four systems. For each system, the program was 
run ten times from one body to ten bodies. One hundred solution points were generated, for each run, with an 
output time-step of 0.01 seconds and a tolerance of 10~^. Total CPU time for solving four systems is calculated 
in terms of CPU time per residual-call or CPU time per Jacobian-call. Total CPU time per residual-call versus 
the number of bodies is plotted in Figure 2. 


i 
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Table 1: Analytical count of operations required for solving three systems with analytical Jacobian 
matrices. 
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7 REMARKS 


The mathematical model of an n-body pendulum with spherical joints and its solution method have been 
presented. The mathematical model derived here is a mixed system of differential and algebraic equations 
(DAE’s) in implicit form. A model of state-space equations can be derived if we do further complex substitu- 
tions. However, in the form of DAE’s the equations of motion provide more sparsity, and this sparsity can be 
exploited when numerical solution methods are applied. The modeling and formulation of an n-body pendulum 
is the first study, for its similarity to a robotic manipulator in structure and for its simplicity in equations. 

We also present four solution methods for solving the equations of motion of this n-body pendulum. The 
first method solves the linear system, J AX = E, directly by LU factorization where J is approximated 
by numerical difference. The second method solves the linear system directly by LU factorization, but J is 
evaluated by its analytical Jacobian matrix with some structural consideration. 

The third method employs the technique of dimension reduction to transform the full system to its reduced 
system, Jr AV = E/*, such that LU factorization is performed on the reduced system only. This technique 
utilizes the sparsity of the system completely and saves a considerable number of computations. The four 
methods also employ the technique of dimension reduction to reduce the size of the system. However, they 
further exploit the structure of the system such that the reduced Jacobian generated possesses a very special 
structure which can be utilized to solve the system in a nearly-linear computational cost. 

The comparative study in computational cost for these four systems in the paper gives strong evidence of 
the success of the theory developed. 


100 



References 


[1] J. C. K. Chou, K. Singhal, and H. K. Kesavan. Multi-body systems with open chains: Graph-theoretic 
model. Mechanism and Machine Theory , 21(3):273-284, 1986. 

[2] W. W. Armstrong. Recursive solution to the equations of motion of an n-link manipulator. Proceeding of 
5th World Congress on Theory of Machines and Mechanisms , 2:1343-1346, July 1979. 

[3] L. R. Petzold. Differential/algebraic equations are not ODE’s. SIAM J. Sci. Stat . Compute 3(3):367-384, 
September 1982. 

[4] L. R. Petzold. A description of DASSL: A differential/algebraic system solver. Scientific Computing , pages 
65-68, North-Holland, 1983. 

[5] Jack C. K. Chou. Modeling, Formulation, and Solution Scheme for an N-Body Pendulum . Technical Report 
No. TR-89008, Program in Engineering Science, Erik Jonsson School of Engineering and Computer Science, 
University of Texas at Dallas, Richardson, Texas, 1989. 

[6] Jack C. K. Chou. Computer-Aided Design Methods for Three-Dimensional Constrained Mechanical Sys- 
tems. Ph.D. Dissertation, Department of Systems Design Engineering, University of Waterloo, Waterloo, 
Ontario, Canada, 1988. 

[7] Jack C. K. Chou. Quaternions, Finite Rotation, and Dynamics . Technical Report No. TR-89007, Program 
in Engineering Science, Erik Jonsson School of Engineering and Computer Science, University of Texas at 
Dallas, Richardson, Texas, 1989. 


101 



